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Abstract 

We deduce the cosmogenic neutrino flux by jointly analysing ultra high energy cosmic ray 
data from HiRes-I and II, AGASA and the Pierre Auger Observatory. We make two determina- 
tions of the neutrino flux by using a model-dependent method and a model-independent method. 
The former is well-known, and involves the use of a power-law injection spectrum. The latter 
is a regularized unfolding procedure. We then use neutrino flux bounds obtained by the RICE 
experiment to constrain the neutrino-nucleon inelastic cross section at energies inaccessible at 
colliders. The cross section bounds obtained using the cosmogenic fluxes derived by unfolding 
are the most model-independent bounds to date. 
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1 Introduction 



Neutrino cross sections at high energies could herald new physics because several extensions of 
the Standard Model predict enhanced cross sections. Examples include electroweak instanton 
processes [T] , black hole production [2] and exchange of towers of Kaluza-Klein gravitons JHj in low 
scale gravity models 0j, and TeV-scale string excitations 5 . The highest center-of-mass energy 
at which the neutrino-nucleon cross section has been measured is about 300 GeV at the HERA 
accelerator. The only practical way to probe this cross section at center-of-mass energies above 100 
TeV is by detecting the interactions of neutrinos with energy above 10 10 GeV incident on Earth. 
By ultra high energy one often means energies above 10 9 GeV, but this usage is not standard. 

While neutrinos with energies above O(10 4 ) GeV have not been observed so far [H], they are 
expected to accompany ultra high energy cosmic rays (that have been observed with energies 
exceeding 10 10 GeV) because almost all potential cosmic ray sources are predicted to produce 
protons, neutrinos and gamma rays with comparable rates. Experiments indicate that the highest 
energy cosmic rays are primarily protons 7 . A guaranteed source of neutrinos, the cosmogenic 
neutrinos, arises from the inelastic interactions of cosmic ray protons on the cosmic microwave 
background [H] , prominently jry — > A + — > mr + , followed by pion decay |U] . The threshold energy 
for this reaction is the Gresein-Zatsepin-Kuzmin (GZK) energy, Eqzk ~ 4x 10 10 GeV 8 . At such 
high energies, the gyroradius of a proton in the galactic magnetic field is larger than the size of 
the Galaxy and it is therefore expected that these cosmic ray protons are of extragalactic origin. 
Since the attenuation caused by the above reaction has a length scale of about 50 Mpc |101 lllj . 
a strong suppression in the cosmic ray spectrum is expected above Eqzk- There is evidence for 
the GZK cutoff in Fly's Eye/HiRes QH H3] data but not in AG AS A [H [To] data. The Pierre 
Auger Observatory |16j (which we will refer to as Auger in what follows) is expected to resolve this 
conflict. 

Estimates of the cosmogenic flux are very model-dependent and differ by about 2 — 3 orders 
of magnitude [111 1171 ITS] . As we discuss below, the uncertainty in the cosmogenic flux is further 
exacerbated by recent evidence from the HiRes experiment that protons dominate the flux at about 
5 x 10 8 GeV jTUj, two orders of magnitude below Eqzk- 

The uncertainty in the cosmogenic neutrino flux must be dealt with before progress can be 
made to extract or constrain the neutrino-nucleon inelastic cross section using cosmic neutrinos. 
We prefer to use cosmic ray data to infer the cosmogenic neutrino flux, thus side-stepping theoretical 
modelling. So long as a lower bound on the cosmogenic flux is established, it will be possible to 
place an upper bound on the neutrino-nucleon inelastic cross section. 

In Refs [201 I21j . estimates of the cosmogenic neutrino flux have been made from separate 
analyses of AGASA and HiRes-II data combined with their predecessor collaborations, Akeno and 
Fly's Eye, respectively. These results have been used to constrain the neutrino-nucleon cross section 
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in Ref. Our approach is different. 

It has been shown in previous work that the apparent disagreement between AGASA and HiRes 
data could be a result of large systematic uncertainties in the energy determinations of the two 
experiments [23]. The energy scale uncertainty of the AGASA experiment is 30% |15j . while that 
of the HiRes detectors is about 20% above 10 10 5 GeV and about 25% above 10 9 5 GeV (2U - The 
data from Auger have an energy uncertainty of about 25% 25 . 

Rather than perform separate analyses of these datasets, we take into account the large un- 
certainties in the energy measurements and perform a combined analysis of Auger j^j, HiRes-I 
(sample collected between June 1997 and May 2005) and II (sample collected between December 
1999 and May 2003) ^3], and AGASA data |14| . We expect the double-counting of events that are 
common to the HiRes-I and HiRes-II datasets to have a negligible effect on our results. We suppose 
that the energy uncertainties for HiRes-I and II are the same and allow a total of 3 floating energy 
scale parameters: one each for Auger, HiRes and AGASA. 

It was thought that the ankle at 10 10 GeV indicates the onset of the dominance of a higher 
energy flux of particles over a lower energy flux, such as an extragalactic component starting to 
dominate over the galactic component of the cosmic ray flux. The dip structure in the vicinity 
of the ankle can be explained as a consequence of pair production from extragalactic protons on 
the CMB |26| . The latter process has a threshold of 10 8 ' 6 GeV thus permitting the interpretation 
that protons dominate even below this energy. Recent evidence for a steepening of the spectrum at 
about the same energy |27| . the second knee, suggests that a low energy transition is consistent with 
the data. Moreover, the HiRes collaboration finds a drastic change of composition across the second 
knee, from about 50% protons just below this knee to about 80% protons above [T2J|. Since the 
end-point of the galactic flux is thought to be comprised of heavy nuclei, a changeover to protons 
is interpreted as the onset of the dominance of the extragalactic flux. We do not analyse events 
with energy below the second knee since they are probably not of cosmological origin. However, 
since composition measurements are strongly model-dependent (via the theoretical uncertainty in 
predicting electron and muon shower sizes and the depth of shower maximum for hadronic showers) , 
the only statement about the transition energy that can be made with some confidence is that it 
lies in the interval 10 8 GeV to 10 10 GeV. 

2 Modelling 

We assume that all observed cosmic ray events above the second knee are due to protons and that 
cosmic ray sources are isotropically distributed. We follow Ref. |20j and write the differential flux 
of particles of type b (b may be protons, p, or cosmogenic neutrinos, j), with energy E arriving 
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at earth as 

d 4 Nh f'°° 

where 

p [l + z(r)] n @(z - z min )@(z max - z) , (2) 

Nb is the number of b particles and A, t and denote area, time and solid angle, respectively. 
This equation needs some explanation. Pb(E; Ei,r) is the probability of a 6 particle arriving at 
earth with energy above E if a proton with energy Ei was emitted from a source at distance r |28| . 
These propagation functions are available from Ref. [23]. Here, po is the comoving number density 
of sources. The step functions restrict the sources to have redshifts between z m i n and z max . We 
take z m i n = 0.012 (corresponding to 50 Mpc) and z m ax = 2. Redshift evolution of the sources is 
parameterized by n, which also mimics changes in z max . Reference [2^ has emphasized that since 
the transition energy is lower than the threshold for pjcMB absorption, sources at high redshift 
are also sampled. Thus, source evolution must be accounted for. For pure redshifting, n = 3. 
Since dz = (1 + z)H(z)dr (where H(z) is the Hubble parameter), Jj, depends on the cosmology 
describing our universe. We adopt a flat cosmological constant-dominated universe with £l m = 0.3 
and H(0) = 71 km/s/Mpc, since these values were used to calculate the propagation functions; 
our results are insensitive to variations in the cosmological model chosen. I(Ei,t) is the injection 
spectrum. 

Throughout, we assume that extragalactic magnetic fields are smaller than 10 -9 G and therefore 
neglect synchroton radiation of protons. 

We employ two methods for obtaining the proton injection spectrum (which we assume to be 
identical for all sources). 

2.1 Power-law injection spectrum 

In the first method, we adopt a power-law spectrum with index a, 

I(E i ,t)=I E?Q(E max -E i ), (4) 

where Iq the normalizes the injection spectrum and E max is the maximum injection energy achiev- 
able through astrophysical processes which we set equal to 3 x 10 12 GeV. The overall normalization 
PqIo is determined from cosmic ray data. 

After extracting a from cosmic ray data, we can determine the cosmogenic neutrino flux re- 
sulting from the GZK mechanism by setting b = v in Eq. (^Q). (We either fix n or determine it 
simultaneously with a). 

For a description of the statistical procedure see Appendix A. 



G b (E,E l ) = ±-J o C °dr 



dP b {E;Ei,; 



dE 
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2.2 Unfolding the injection spectrum 

The essential idea behind deconvolving the injection spectrum from cosmic ray data is described 
in Ref. [H]. 

The proton injection spectrum (with statistical uncertainties) is obtained by inverting the ob- 
served proton spectrum. Use must be made of the propagation function of the proton. Although 
this function is not invertible in general, the injection spectrum may be unfolded under the as- 
sumptions made earlier that cosmic sources are isotropically distributed (within a redshift range, 
[z m in, z max \) and the redshift evolution can be parameterized by (1 + z) n . Then, Eq. (Q) takes the 
form of a matrix equation after writing the integral as a sum over energy, 

I = G^3 b , (5) 

where I is the discrete version of the injection spectrum I(Ei,t). Since cosmic ray data give Jb 
with b = p, we can find I. This seemingly simple procedure is fraught with technical difficulties. 
We relegate the details of the statistical methodology to Appendix B. 

It is noteworthy that this method is applicable even for AGASA data that do not show evidence 
for a GZK suppression; events beyond the GZK suppression can not be accounted for by a single 
power-law injection spectrum. 

The cosmogenic neutrino spectrum (with statistical uncertainties) is inferred from the proton 
spectrum by using 

J„ = G„I. (6) 

A lower bound on the neutrino flux is obtained by assuming that only the events at and above 
the pile-up below the GZK energy are protons. This is consistent with the GZK interpretation that 
the energy of protons close to the GZK energy degrades significantly causing the pile-up. The upper 
bound is found by supposing that all the observed events above the second knee are protons. The 
events above the GZK energy could arise simply because the injection spectrum was sufficiently 
large. 

3 The cosmic ray spectrum and the cosmogenic neutrino flux 

We first consider the case of a power-law injection spectrum. In our analysis we only consider data 
in the energy range 10 9 - 6 GeV to 10 11 GeV. Including data outside this range gives results which 
have a goodness of fit (gof) below O(10~ 10 ). That this is the case for data above 10 11 GeV is 
understood as a consequence of the inability of a power-law injection spectrum to explain super- 
GZK events. The AGASA data point at 10 955 eV is statistically significant and discrepant with 
data from other experiments. 

In Fig. ^ we plot J p vs. E for three different values of n and for the case in which n is a free 
parameter in the fit that is allowed to vary between and 6. The solid line is the best-fit to the data 
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in each case. In all cases, a spectral index close to —2.4 is favored. The shaded band and quoted 
uncertainties correspond to models that are consistent with the data at the 2a C. L. The error bars 
on the data points are la uncertainties. The insets show J p E 3 vs E to enable comparison with 
results presented in other papers. The y-axis units on the left-hand side are GeV _1 cm _2 s _1 sr _1 
and those on the right-hand side are eV m — 2 s sr _1 . Note that since the energy scales of the 
experiments have large uncertainties, it is misleading to plot J P E 3 . From panel (d) it can be seen 
that the 2a uncertainty in n spans almost the entire range within which n was permitted to vary. 
Although n < 1 is disfavored at 2a, we will show results for n = because its exclusion is not very 
significant. 

The goodness of fit of these joint analyses is poor. This is primarily because the initial Auger 
data are noisy; see Tables 1 and 2. For now, we proceed to determine the corresponding cosmogenic 
neutrino fluxes. 





x 2 


dof 


gof 


All data 


86.3 


45 


2.1 x 10" 4 


—Auger 


37.6 


32 


0.23 


-AGASA 


58.8 


32 


2.7 x 1(T 3 


-HiRes 


65.1 


24 


i.2 x i(r 5 



Table 1: The x 2 f° r the number of degrees of freedom (dof) and the corresponding goodness of 
fit (gof) for analyses with n = 3. The first row is for an analysis of all data. The cosmic ray 
spectrum is shown in the panel (b) of Fig. ^ Each subsequent row is the result with the indicated 
dataset removed from the analysis. Removal of the Auger data from the analysis improves the gof 
considerably. 





x 2 


dof 


gof 


Auger 


38.7 


11 


6 x 10" 5 


AGASA 


10.4 


11 


0.49 


HiRes I+II 


16.5 


19 


0.62 



Table 2: Similar to Table 1, except that the results are for separate analyses of each experiment. 
The analysis of the Auger dataset yields a poor gof. 

With the cosmic ray spectrum in hand we can compute the neutrino flux J v produced in the 
GZK chain by using Eq. Q with b = v. The resulting cosmogenic neutrino fluxes (summed 
over flavors) are shown in Fig. [21 We do not show the case in which n is free, since it fills the 
region between the n = and n = 6 bands. Again, the bands correspond to the 2a C. L. The 
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95% C. L. model-independent upper bound obtained by the RICE collaboration from their latest 
data compilation is also shown. 

We have seen that a combined analysis including Auger data has an unsatisfactory goodness 
of fit even if super-GZK events seen by AGASA are excluded. This suggests that an alternative 
method to obtain the injection spectrum is called for. 

The injection spectra obtained by unfolding are shown in Fig. |21 The widths of the bands are 
the 2a C. L. determinations. There are four different spectra since the unfolding was carried out 
experiment-by-experiment; we have not attempted an inversion of the joint data because there is 
too much freedom, rendering the results meaningless. So that a straightforward comparison can 
be made with the power-law injection case, we have analysed data in the energy range 10 9 ' 6 GeV 
to 10 11 GeV; see panels (a) and (b) of Fig. |21 However, since the unfolding method is applicable 
even for super-GZK events, we have performed a separate analysis in which we include data up to 
3 x 10 11 GeV; see panels (c) and (d) of Fig. 

In the unfolding procedure it is assumed that there is no flux beyond the energy range spanned 
by the data. The sharp edges at the two ends of the spectra result because the spectra can not 
be extrapolated outside the energy interval. This is not true for the case in which a power-law 
injection spectrum is assumed. The bold line corresponds to a power law injection spectrum with 
spectral index a = —2.5 in the n = plots and a = —2.3 in the n = 6 plots. The value of a is 
pinned-down by the statistically significant data at the lower end of the energy range. Thus, we 
have used the same spectral index in panels (a) and (c) and in panels (b) and (d). Note that the 
bold lines are extrapolated outside the range covered by the unfolding method. Results for values 
of n between and 6 are intermediate to those in the left-hand and right-hand panels of Fig. |3J 

The neutrino spectra in Fig. 0] are obtained from Eq. The bands correspond to the 2a 
C .L. The bold solid lines in the top and bottom panels are for power-law injection spectra with 
E m ax = 10 11 GeV and E max = 3 x 10 11 GeV, respectively. The bold dashed lines are the neutrino 
spectra of Fig. [2] and are provided for comparison. There we had set E max = 3 x 10 12 GeV. The 
agreement between the two methods is striking and provides support for the assumption that the 
injection spectrum has the form of a power-law. The poor fit obtained in the combined analysis 
with a power-law spectrum is likely due to systematic uncertainties that have not been accounted 
for. 

Since the flux obtained from the unfolding procedure can be sizeable, it is essential to check 

that the accompanying cosmogenic photon flux is not in conflict with the EGRET observation of 

the diffuse gamma ray flux [3U] , as emphasized in Ref . [3^ . We evaluated the photon flux using the 

publically available software, CRPropa [32| . and find comfortable consistency with the EGRET 

bound. This is because the injection spectra we have found fall steeply with energy and because 

1 Incidentally, we first confirmed that the cosmic ray spectra generated using the propagation functions |29| and 
CRPropa are identical. 
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the largest value of E max we consider is 3 x 10 11 GeV. As a consequence, the contribution to the 
photon flux measured by EGRET is tiny. 

4 Constraints on the neutrino-nucleon cross section 

In the energy range of interest (10 8 to 10 11 GeV) the strongest bounds on the neutrino flux are 
those of the RICE experiment [33] • (The current bound from the ANITA experiment is stronger 
than that from RICE only above 10 11 GeV |34j ) . The RICE bound is still much weaker than the 
determinations we made in the previous section. This allows us to place an upper bound on the 
neutrino-nucleon cross section; new physics can not increase the Standard Model cross section too 
much or high energy neutrinos would have been observed at neutrino telescopes. Note that if the 
upper limit of the cosmogenic neutrino flux from either method had been above the RICE limit, 
we could also have placed a constraint on the minimum required suppression of the cross section 
at ultra high energies. 

The RICE collaboration has employed a method to obtain neutrino flux limits that are inde- 
pendent of specific neutrino flux models. See Appendix II of Ref. }331 . So long as new physics does 
not alter the energy dependence of the Standard Model cross section drastically (and only changes 
the overall normalization), we can constrain the neutrino-nucleon inelastic cross section by simply 
taking the ratio of the RICE bound to that of the lower bound of the cosmogenic neutrino flux. 

Our 95% C.L. upper bound on the cross section (in units of the Standard Model cross section) 
is shown in Fig. [5] for n = 0, 3 and 6. In each case, the colored curves correpond to the limit 
obtained by using the propagation inversion procedure and the black curve corresponds to a power- 
law injection spectrum. Bounds derived from RICE are valid only for a v < 1 mb Our cross 
section bounds are applicable only in the unshaded region. 

5 Summary 

We obtained the proton injection spectrum from cosmic ray data using two methods. In the first, 
we performed a combined analysis of HiRes-I and II, AGASA and Auger data in the energy interval, 
10 9 ' 6 GeV to 10 11 GeV, using a power law injection spectrum. We found the quality of the fit to 
be unsatisfactory even though we did not include super-GZK events in the analysis. The primary 
reason for the poor fit is apparently that the initial Auger data are noisy. In the second method, 
we implemented a regularized unfolding of the injection spectrum for each dataset separately. 

We found the resultant cosmogenic neutrino flux using the injection spectra from both methods 
to be in excellent agreement with eachother. This suggests that the injection spectrum is well- 
modeled as a power-law and that the poor fit mentioned above is due to unaccounted-for systematic 
uncertainties. 



8 



Using the model- independent limit on the neutrino flux from RICE, we constrained the neutrino- 
nucleon cross section under the assumption that new physics modifies the charged-current and 
neutral-current interactions by the same constant factor. 

Our results are succintly summarized in Figs. 
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A x 2 analysis 

The observed flux J± is given by rap/ej, where nf and ej is the binned event rate and effective 
exposure, respectively, in bin i of an experiment. The theoretically predicted flux J\ is obtained 
from Eqs. (|H4j) . It depends on the overall normalization pqIq, the spectral index a and on redshift 
evolution as parameterized by n. The theoretical event rate n\ is related to the theoretical flux by 



In order to incorporate the effect of the energy scale uncertainty on the event rates, we replace 
the (observed) energy Ei by (1 + j)E{, where 7 is the fractional energy scale uncertainty of the 
experiment. The details follow the implementation in GLoBES [33]. Thus, the theoretical flux 
becomes a function of polo, ot, n and 3 7's. 

Since the the observed event rate can be very low or even zero, we use the Poissonian x 2 -function 



It is crucial to propagate the errors and their correlations consistently as we are interested in both, 
the proton flux and the derived cosmogenic neutrino flux. The supposition that the errors are 
uncor related would lead to an incorrect result. The allowed region at confidence level CL is defined 
by requiring that A% 2 < XcL' w here Xcl * s gi yen by the CL th percentile of the x 2 -distribution 




(see, e.g. jSH]) 
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after accounting for the 5 or 6 free parameters in the analysis 2 . We use a Markov Chain Monte 
Carlo (MCMC) method [3ZJ based on the Metropolis-Hastings algorithm. Specific care has to be 
taken to ensure that the chain has reached equilibrium, i.e., its asymptotic state. We employ 
the convergence diagnostic of Ref. |38| since it is conceptionally simple and does not incur a large 
computational burden. 

B Unfolding procedure 

The propagation of ultra high energy cosmic rays is described by Eq. ((J) which is a Fredholm 
equation of the first kind. The authors of Ref. j2U] have attempted to infer I(Ei) from the observed 
flux Jb(E) and Gb(E,Ei) via direct inversion. Unfortunately, the problem is ill posed unless regu- 
larized unfolding is carried out; see Ref. [35] for a detailed exposition. An illustration of this fact 
is provided in Fig. H3 

We follow Tikhonov regularization described in Ref. [3U|. The x 2 definition of Eq. Q must be 
modified to include a regularization function <S(n*) and a regularization parameter [3, 

xV)+/^(n'), (8) 

with the understanding that the n\ are the parameters in the fit; i.e., the n\ are not derived from 
a model. It can be shown that ^-minimization with (3 = is an unbiased estimator with the 
smallest variance but which yields highly oscillatory solutions. One the other hand, the (3 — > oo 
case yields a maximally smooth estimator with vanishing variance and a clear bias (since the result 
no longer depends on the data) . The parameter (3 determines the relative weight placed on the data 
in comparison to the degree of smoothness of the solution. In a Bayesian spirit this is equivalent 
to assuming a prior which favors smooth solutions with (3 controlling the width of the prior. The 
smaller (3 is, the less impact the prior will have and vice versa. 

It remains to choose <S'(n*) (to obtain a smooth solution) and (3 (to define the trade-off between 
bias and variance). For our application, the mean square of the second derivative 

N-2 

5(n*) = -J2 {-< + 2 <+i " < + 2? , (9) 
i=l 

works very well. 

Our criteria for obtaining (3 are that injection spectrum be positive definite and that the overall 
X 2 for an individual dataset not increase by more than 50% relative to that obtained from the 
power-law fit. Values of (3 above unity satisfy the former requirement, while (3 needs to be less 
than 10 to satisfy the latter. We select [3 = 5 for all four datasets, which adequately suppresses 
any spurious oscillations. 

2 3 energy scales, 1 normalization, 1 spectral index and optionally the evolution parameter n. 
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We estimate the variance of the result by generating 1000 random realizations of the data sets 
by assuming a Poisson distribution in each bin with a mean value nf. For each of these random 
realizations we repeat the above procedure and construct the covariance matrix C which is then 
used to obtain upper and lower bounds on the proton and neutrino fluxes by retaining only those 
realizations which satisfy 

n*C-V) T <X&L- (10) 
Here, Xcl 1S gi ven by the CL th percentile of the values obtained for the left hand side. 
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Figure 1: The flux of J p of ultra high energy cosmic rays that simultaneously fits Auger, HiRes-I 
and II, and AGASA data in the energy range 10 96 GeV to 10 11 GeV assuming a power-law injection 
spectrum with spectral index a. Redshift evolution of the sources is parameterized by a power-law, 
(1 + z) n . The shaded bands and the quoted uncertainties are at the 2a C. L. and the error bars on 
the data points are la uncertainties. The energy scale of each experiment has been allowed to vary 
within the experiment's energy scale uncertainty. The insets show J p E 3 to facilitate comparison 
with other analyses. 
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Figure 2: The 2<t C. L. determinations of the cosmogenic neutrino flux (summed over all flavors) 
corresponding to the injection spectra of Fig. ^ The light shaded region is the 95% C. L. model- 
independent upper bound from RICE. 
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Figure 3: The 2a unfolded injection spectra from Auger, HiRes-I and II, and AGASA data. The 
upper (lower) two panels are for an analysis of data with energy between 10 96 GeV to 10 11 GeV 
(10 9 - 6 GeV to 3 x 10 11 GeV). Redshift evolution of the sources is parameterized by a power-law, 
(1 + z) n . The bold line in the n = (n = 6) plots is the injection spectrum of the form E^ 2 ' 5 



17 




Figure 4: The 2a C. L. determinations of the total cosmogenic neutrino flux corresponding to the 
injection spectra of Fig. |21 The bold solid lines in the upper and lower panels are for power-law 
injection spectra with E max = 10 11 GeV and E max = 3 x 10 11 GeV, respectively. The agreement 
between the two methods is very good. The bold dashed lines are the neutrino spectra of Fig. O 
which were obtained for power-law injection spectra with E max = 3 x 10 12 GeV. 
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Figure 5: 95% C.L. upper bounds on the neutrino-nucleon cross section using our determinations 
of the cosmogenic neutrino flux and the model-independent flux limit from RICE. Our bounds are 
only valid in the unshaded region for which a v < 1 mb. 
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Figure 6: Two different injection spectra, I\{Ei) and ^(-E/j), that yield the same cosmic ray spec- 
trum J p (E) on folding in G p (E, E^). A naive unfolding of data corresponding to J P (E) (i.e., with 
random fluctuations included) could result in a highly oscillatory injection spectrum like ^(-^i) 
which may be not be positive definite at all energies. The goal of regularized unfolding is to extract 
the smooth spectrum I\(Ei) with the minimal introduction of bias. 
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